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ABSTRACT 

In this paper we report on the formation of magnetically-levitating accretion disks around supermas- 
sive black holes. The structure of these disks is calculated by numerically modelling tidal disruption 
of magnetized interstellar gas clouds. We find that the resulting disks are entirely supported by the 
pressure of the magnetic fields against the component of gravitational force directed perpendicular to 
the disks. The magnetic field shows ordered large-scale geometry that remains stable for the dura- 
tion of our numerical experiments extending over 10% of the disk lifetime. Strong magnetic pressure 
allows high accretion rate and inhibits disk fragmentation. This in combination with the repeated 
feeding of magnetized molecular clouds to a supermassive black hole yields a possible solution to the 
long-standing puzzle of black hole growth in the centres of galaxies. 



1. INTRODUCTION 

It is believed the growth of supermassive black holes 
(SMBHs) in centres of galaxie s is enabled by gas accre- 
tion from surrounding disks ( |Lynden- Bell 1969) which 
have been observed with increasing precision by mo dern 
telescopes ( jMiyoshi et~aT1|1995l [Jaffe et al . j | 2QQ4|) . In 
the early theoret ical work ( |Lynaen-B ell 1969; Sha kura| 
' fc Sunyaevj|1973| ) it has been suggested that that mag- 
netic stresses play an important role in driving the accre- 
tion by enabling the outward angular-momentum trans- 
port through the disk. This s uggestion has been put on 
a firm theoretical footing by Bal bus fc Hawley] ( |1991[ ) 
discovery of the importance of magnetorotational insta- 
bility (MRI) in astrophysical disks, and by the subse- 
quent work, that demonstrated the ability of MRI to 
build an d maintain substantial magnetic stresses inside 
the disk ( Brande nburg et all 19951 |Stone et al."|l996[ IhT 



rose et al.||2(J(J6| |Davis et al.||2(J!U| ) . All of the numerical 
studies to date have demonstrated MRI-generated mag- 
netic stresses which are associated with the sub-thermal 
magnetic fields in the disk mid-plane. 

One of the central unresolved issues of feeding SMBHs 
has been the tendency of all modelled ex tended gaseous 
disks to clump due to their self-gravity (|Kolykhalov fc| 
Syunyaev 1980; S hlosman fc Begelman||19871 |Shlosman[ 
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2009). Such chok- 



mg ol the accretion flow is a maj or obstacle in SMBH 
growth. It has been conjectured jShibata et al. 1990; 



Mac hida et"aLl|2000l |Pariey e t al. 2003; Mac hida et al. 



2006 ; |Begelman fc Pringle||2007| |6da et al.||2009[) that 



m some astrophysical disks magnetic stresses may be- 
come dominant relative to the mid-plane gas pressure, 
and that these disks may effectively resist fragmentation. 
In this paper we investigate the formation of accretion 
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disks by performing numerical simulations of collisions 
between magnetized gas clouds and a black hole. It has 
been suggested that such collisions may be responsible 
for feeding the supermassive black holes at the centers 
of galaxies ( King fc Pringle|2007 Wardle & Yusef-Zadeh 



2008 2012[ ) and that it may have lead to the formation' 
ol the stellar disc in our own Ga l actic Center (|Sanders| 
T9981 ILevin fc Beloborodovl[2003l fPaumard et aT] f2006 ; 
Bonnell fc Rice|2008|IHobbs fc Nayakshin|2009p . We find 
that the resulting disks are completely dominated by the 
magnetic field pressure, and display high accretion rates 
due to the Maxwell stress associated with the large-scale 
magnetic field the structure of which remains stable over 
the duration of the simulation. The Toomre-Q factors of 
these naturally-formed magnetically-levitating accretion 
disks (MLAD) indicate their stability to gravitational 
fragmentation. Therefore, MLADs represent a new class 
of accretion-disk solutions which may play an important 
role in feeding the supermassive black holes. 

2. SIMULATIONS SETUP 

We model a collision between a magnetized gas cloud 
and a SMBH using a new moving-mesh ideal MHD 
scheme (see Appendix). We choose an equation of state 
^Pgas = CgP, where c s = 0.03 vk] here vk is Keplerian ve- 
locity around a SMBH. The temperature in this setup is 
T = 1.63 x 10 4 K(0.1pc/i?), which, in absence of mag- 
netic fields, would produce disks with Hq/R = 0.03. In 
the presence of magnetic fields the effective scale-height 
is modified by the magnetic pressure, H = Ho ^/l + 
where f3~ x = P m /P g is the ratio of magnetic to gas pres- 
sures. In order to isolate the effects of magnetic fields on 
the disk formation process, we ignore effects of the gas 
self- gravity in our calculations. 

We conduct three high resolution simulations with 
1.6 x 10 7 particles in which a magnetized gas cloud with 
mass and radius 3.5 pc and 8.8 x 1O 4 M , respectively, 
is collided with a 3.5 x 10 6 M S MBH. The geomet ry of 
the initial conditions is taken from Alig et al. (2011 ) and 
repeated in Table [I] in particular, we use initial condi- 
tions from their simulations COl, VOl and V02, where 
a molecular cloud has an impact parameter of 2 pc, 3 pc 
and 3pc respectively. On top of these initial conditions 
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Fig. 1. — Snapshots of density structure in the XY plane at different times for V01 model. The times are shown in units of 0.047 
million years, which corresponds to 0.0, 0.047, 0.096 and 0.240 million years for top-left, top-right, bottom-left and bottom-right panels 
respectively. The unit of length is a parsec, and unit of density is 4.1 • 10 6 cm -3 m u . 



TABLE 1 

This table shows geometry of initial conditions: Column 1 
(id) is the name of a run, column 2 (r c \) show adius of a 

cloud, Column (v x ) is cloud infall speed in km/s, and 
finally last coumn (6) is the cloud 's impact parameter. 



ID 


Rd [pc] 


v x [km/s] 


b [pc] 


C01 


3.5 


120 


2 


V01 


3.5 


30 


3 


V02 


3.5 


50 


3 



we also impose a uniform magnetic field that threads 
both the cloud and vacuum regions. The magnetic field 
strength is such that the resulting magnetization inside 
the cloud is f3 = 1, which corresponds to \B\ « 100 /iG 
and dimensionless mass-to-flux ratio C ~ 5, where C = 
(M/^>)/ A /5/(97T 2 G) (|Mouschovias fc Spitzer||l976| |Mac| 



Low fc Klessen|20'04j ). This tield strength c orresponds to 
the large-scale held m the Galactic Centre flYu sef-Zadeh 
Mprris|1987||Morris Yusef-Zadeh|1989||Crocker et a 



2010p . The initial magnetic field orientation was such 
that each of the components of magnetic field have the 
same magnitude, namely B x = B y = B z = B/y/3. The 



vacuum is modelled with fluid 10 6 times less dense than 
the cloud density (n = 0.01cm -3 ), which we also use it 
as a floor density to avoid local density contrasts larger 
than ~ 10 T that our code cannot deal with due to use of 
single precision floating point arithmetics. 

The computational domain is a periodic box with 
32 x 32 x 32 pc 3 volume, which is large enough not to 
influence physical processes occurring in sub-parsec re- 
gions. The mass and distance units were [M] = 10 5 M 
and [R] = 1 pc respectively, which sets the time unit 
[T] w 0.047 Myr, magnetic field units [B] « 5.40 mG, 
and the speed unit [V] ~ 20 km/s. The simulation lasted 
till T end = 5.0, which corresponds to 0.24 Myr or 4.7 
orbital periods at R = lpc. The inner boundary con- 
ditions are applied only within 0.02 pc from the SMBH, 
which we regard as the inner disk boundary, as follows. 
Any particle within 0.01 pc is removed from the compu- 
tational domain, and in the transition region between 
0.01 pc and 0.02 pc we set the density to the floor value, 
and both the velocity and the magnetic field to zero. 



3. RESULTS 
3.1. Geometry of the collision 





Fig. 2. — Density snapshot in the XY plane at the of V02 (top left-panel) and C01 (top right-panel) simulations. As in Fig.^ the unit 
of length is a parsec, and density is shown in units of 4.1 • 10 6 m u /cm 3 . The bottom left and right panels show density profile in XY plane 
at the end of V01 simulation and the associated divergence error, divB = | V • B|/i/|B| where h is cell size, respectively. 

A collision between a molecular cloud and a supermas- 
sive black hole is a violent event occurring on dynami- 
cal time-scale. Since fluid elements generally have non- 
zero angular momentum, the natural outcome of such an 
event is a formation of a disk. Hydrodynamical simula- 
tions of such collision event robustly show a formation 
of an eccentric disk around SMBH with the d isk geome 



try b eing dependent on the initial conditions (Alig et al. 
2Q11[). Our aim in this work is to study similar event 



but in strongly magnetised regime in which the initial 
magnetic pressure in the cloud is in equipartition with 
the gas thermal pressure. 

In Fig. [I] we show snapshots of the gas density in the 
XY plane at t=0, 47, 94, and 240 thousand years. In 
the first hundred thousand years, the cloud experience a 
violent collision with the black hole. In particular, the 
bow-shock, which can be seen in the bottom-left panel 
as a large curved region with density jump just above 
the disk, is formed by isothermal shock guards the newly 
formed inner disk from the destructive effect of the in- 
coming fluid. The outcome of this collision event is a 
formation of a parsec-size gas disk with irregular density 
structure. Similar disks where formed in other simula- 
tions as can be seen in the top two panels of Fig. [2] This 



can be contrasted with Alig et al. ( 2Q11[ ) where simula- 
tions C01 and V01 have hnal differently shaped gas disks. 
Finally, in the bottom right panel of Fig. [2] we show the 
divergence error, |V • B|/i/|B where h is the cell size, in 
the final snapshot of V01 simulation. 

In Figj3|we show magnetic field geometry in different 
regions of the disk at the end of the V01 simulation, 
which corresponds to approximately 150 orbital periods 
at R = 0.1 pqj The top- left panel show magnetic field 
lines originated in the central region of the disk and ex- 
tend above and below mid-plane. The magnetic field in 
this regions is dominated by poloidal components. The 
top right panel shows mid-plane magnetic field structure 
in the central region (R ~ O.lpc). The field lines ap- 
pear regular and tightly winding in azimuthal direction, 
which is result of strong Keplerian shear inside the disk. 
In the bottom-left panel, we also show magnetic field in 
the mid-plane region but further away from the center 
(R ~ 0.5 pc). While magnetic field is still stretched in 
azimuthal direction, in contrast to the central regions it 
shows less regular structure. In the bottom-right panel 

6 Since the inner region of the disk is formed at approximately 
t ~ 0.05 million years, the actual number of disk revolutions at 
R ~ O.lpc is <100. 
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Fig. 3. — This figure shows magnetic field geometry in different regions of the disk. The density is shown in the XZ plane, with values 
given in log(p/(4.1 • 10 6 cm 3 /ra u ). The top two panels show magnetic field geometry in the central region of the disk (R ~ 0.1 pc), whereas 
two bottom panel focus on outer regions of the disk (R ~ 0.5 pc). While magnetic field geometry is mostly uniform, the gas density show 
irregular structure. 



we show magnetic field in the disk corona, where mag- 
netic field shows regular large-scale azimuthal pattern. 

3.2. Vertical structure 

This section focuses on the vertical disk structure in 
our simulations. In particular, we studied vertical struc- 
ture of the disk at R w 0.1 pc, which is far enough from 
the inner boundary and close enough that the disk per- 
formed approximately 100 orbital periods by the end of 
the simulation. One of the crucial properties in mag- 
netised disk simulations is the MRI quality factor, Q, 
which is related to the number of resolution points, 
e.g. grid-points, particles or mesh-cells, per MRI fastest 
growing mode wave-length. With this number being 
too low (Q < 8), the simulation may fail to faithfully 
model magneto-rotational instability (e.g. Hawley et al. 
( |2011| )). Due to the nature of our simulations, it was im- 
possible a priori to identify which of our simulations can 
faithfully model long-term disk evolution. As a result, we 
computed MRI quality factors in vertical and azimuthal 
directions at the end of our simulations, and check which 
of the simulations were able to resolve MRI. Namely, we 
compute Q z = X z /h and = A^/Zi, where h is the size 
of a resolution element and X z $ « 2tt\B z ^| / \yfiXl) . 

In Table|2] we show vertically averaged quality factors 
at R ~ 0.1 pc. This table demonstrates that all simu- 
lations have > 8, which means they can faithfully 



TABLE 2 

MRI QUALITY FACTOR FOR EACH OF OUR SIMULATIONS 



ID 


Q z 


Q<\> 


C01 


3 


26 


V01 


11 


60 


V02 


6 


48 



model non-axisymmetric MRI. However, only V01 simu- 
lation qualifies when it comes to axisymmetric MRI, and 
therefore we focus our study of vertical structure on this 
simulation. 

We compute scale- height, H, at this radius by fitting 
an isothermal density profiles in approximately two scale- 
heights. The resulting scale- height is H « 0.01, which 
gives H/R ^0.1 (top- left panel in Fig.|4|. The radial 
temperature dependence is expected to produce disks 
with the scale height H/R = 0.03 y/l + ft' 1 which, for 
p- 1 = Pm/Pg ~ 10 found at R = 0.1 pc, gives 

We studied vertical structure of the disk in V01 sim- 
ulation at R « 0.1 pc, which is far enough from the in- 
ner boundary and close enough that the disk performed 
approximately 100 orbital periods by the end of the sim- 
ulation. We compute scale-height, H, at this radius by 
fitting an isothermal density profile in approximately two 





Fig. 4. — This figure shows vertical structure of the disk model at the end of the best resolve simulation (V01). The top left panel shows 
density (red line with open circle) and the fit of the expected density structure (blue line). The top right panel shows the deviation of 
the azimuthal velocity from the Kerplerian value, the bottom left panel shows Maxwell stress, and the bottom right panel show vertical 
dependence of the ratio of the gas pressure to the magnetic pressure. 



scale- heights. The resulting scale- height is H w 0.01, 
which gives H/R « 0.1 (top-left panel in Fig.[|). The 
radial temperature dependence is expected to produce 
disks with the scale height H/R = 0.03 yjl + ft' 1 which, 
for p- 1 = Pm/Pg ~ 10 found at R = 0.1 pc, gives 
H/R ^0.1, consistent with the simulation data (bottom- 
right panel in Fig.[i|. We also studied the deviation 
of azimuthal velocity form the Keplerian velocity at 
R = 0.1 pc as a function of height, and found that az- 
imuthal velocity variations are less than a percent for 
\z\ < H (top-right panel in Fig.[4|. It is therefore jus- 
tifiable to assume that the disk angular velocity is con- 
stant on cylinders. Finally, in the bottom left panel of 
the Fig. [4] we show Maxwell stress «m = —{BrB^/P^ 
which is approximately 0.1 within the scale- height. 

In Fig. [5] we show azimuthally averaged magnetic field. 
This figure show that the magnetic field is confined 
within few scale- heights of the mid-plane. The magnetic 
field is dominated by the azimuthal component that is 
an order of magnitude larger than the radial one. Ver- 
tical component, B Z1 is much smaller compared to both 



B r and B^ for \z\ < H/2 (in this figure both B r and B z 
strength are magnified by a factor of 10). 

All of our simulations show similar vertical confine- 
ment of the field, which can be interpreted as a result of 
the disk formation: a combination of isothermal shocks 
that amplify magnetic field and Keplerian shear which 
generates strong azimuthal field component. However, 
we would like to stress that the field confinement in our 
best resolved model (V01 ) is in a good agreement with 
Johansen fc Levin] (2008) - hereafter referred to as JL, 
who find similar results in their shearing box models. In 
the JL shearing-box simulations, which were performed 
with a grid-based Pencil Code, the disk field was ini- 
tially in equipartition with the gas pressure, but evolved 
by Parker and magnetorotational instabilities to a mag- 
netic field configuration in the vertical direction similar 
to what we see in our disk which is formed via a colli- 
sion of a magnetized gas cloud with the black hole. It is 
significant that the two simulations that are so different 
in their approach give vertical structure of B^ that is 
in a good agreement with each other. In particular, the 
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Fig. 5. — Vertical dependence of azimuthal (red lines with open 
circles), radial (blue lines with open squared) and vertical (green 
line with filled squares) magnetic field components. The dashed 
and solid lines show low- and high-resolution simulations respec- 
tively. Both radial and vertical magnetic field components are 10 X 
magnified. 




Fig. 6. — Sketch of the evolution of a magnetic field line subject 
to Parker instability, Coriolis force and Keplerian shear. The green 
straight tube in all panels shows the original field line along cp di- 
rection in the mid-plane, the white-blue 3D tube depicts the actual 
field line, and its projection onto the r — cp plane is shown in black. 
The top-left panel shows the undulant distortion of the original 
field line due to Parker instability in the <p — z plane. The fluid ele- 
ments that slide along the field lines (red and blue arrows) towards 
the mid-plane are acted upon by the Coriolis force transforming 
the field line into a helical form. Projection of the field line onto 
r — <p plane shows creation of a radial component (bottom-left). 
Stretching of the radial component by Keplerian shear regenerates 
the azimuthal component (top-right). Finally, subsequent stretch- 
ing of the line generates oppositely oriented magnetic field above 
mid-plane and increases the strength of the azimuthal component 
of magnetic field in the mid-plane (bottom-right). 

azimuthal component changes sign at 2-3 scale-heights 
above the mid-plan^J and the strength of the mid-plane 
field is ten times the value of the reversed field. 

It is likely that the process responsible for flux con- 
finement is of a similar nature as described in JL. We 
sketch it in Fig.[6j The strong azimuthal magnetic field 
is subject to Parker instability. As the fluid elements 
slide along the field line towards the mid-pl ane (top- left) 
they are acted upon by the Coriolis force ( jHanasz et al.| 
2002), and due to this the line becomes helical (bottom- 
left). This generates radial field, which via Keplerian 
shear regenerates the azimuthal field that was lost to 

7 In Fig. 7 of JL the scale-height can be increased by y/2 due to 
extra support provided by the magnetic pressure. 



Parker instability (top-right). Further shearing generates 
oppositely oriented azimuthal field above the mid-plane 
which reconnects with the field of the original orienta- 
tion. This decreases its magnitude or even reverses its 
direction, whereas the field at the mid-plane is amplified 
while maintaining its original direction (bottom-right). 

3.3. Radial structure 

The radial evolution of the disk is driven by the flow 
of matter from the outer to the inner regions. This flow 
is enabled by the effective viscosity from magnetohydro- 
dynamical stresses. Within the viscous time-scale, the 
disk reaches steady radial structure which is computed 
by applyin g conservation laws and theory of steady thin 
disks (e.g. |Frank et aL] [2002). Here we present some 
numerical evidence lor an extra constraint, the conser- 
vation of azimuthal magnetic flux, $ = J dS where 
dS = dz v Y dt. The set of equations that describe disk's 
radial structure is 



M = 2irRY>v v , 
$ = 2HB U 



'if »n 



M = 3n a acc E H 2 ft. 



(1) 
(2) 
(3) 



Here, Eq.j2| describes the frozen-in condition of mag- 
netic fielq^P In what follows we assume that the right 
hand side of these equations are constants for steady- 
state disks. We also assume that accretion viscosity a acc 
is set by Maxwell stresses 



(B r Bp 
4ttP 



(4) 



where P = ft 2 H £ is the total pressure. The disk scale- 
height is 

H y/l 



H 



(5) 



where Ho = c s /Q is a hydrodynamical scale-height. Ac- 
cording to the numerical results from the previous sec- 
tion, the magnetic pressure is entirely dominated by B^, 
which gives P m ~ B^/Stt. 

Using these equations and noting that B^/8tt = 

^ 2 i7S, we derive the radial dependence of disk scale- 
height 



H 
R 



3tt 2 GM(1 + /3) 
40 a m C 2 (nR) 3 



(6) 



where we write 6 = v / 5/(9tt 2 G)M/C. In the limit, f3 <C 
1 the radial dependence of disk magnetization is 



fr 



3ir 2 GM 
40a m C 2 (fti?) 3 



#0 



(7) 



In a general case, Hq is self-consistently computed by 
solving radiative transfer equation in the vertical direc- 
tion. Therefore, the magnetization depends on the ther- 
mal properties of the disk. However, in our simulations 



Hq = 0.03 R from which we have f3 oc 



-2/5 



8 The condition can also be derived from the induction equa- 
tion by consideration of the radial advective flux of the vertically 
integrated azimuthal magnetic field, f dz. 
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These equations are consistent with our simulations 
throughout most parts of the disk (bottom panel in 
Fig.[7|, except in the regions close to the disk inner 
boundary. We also do not expect the model to hold 
for R > 0.4 pc, where the viscous time estimated us- 
ing steady thin disk approximation is much larger than 
the duration of simulation. Agreement with the analyti- 
cal model beyond this radius implies that disk evolution 
there occurs at higher than viscous rate derived from the 
steady thin disk theory. We also note that the gas density 
distribution in the disk is not steady, but exhibits clumpy 
and filamentary structures (right panel in Fig.[3|. This 
is also reflected in the irregularity of the surface density 
profile in Fig.[7| Nevertheless, agreement of radial depen- 
dence between the model and simulations indicates that 
azimuthal magnetic flux is conserved in MLADs during 
accretion. 



TABLE 3 

Accretion rate and effective viscosity 



0.001 



R[pc] 


M[M /yr] 


«acc 


a 


0.05 


0.035 


0.06 


0.03 


0.1 


0.066 


0.09 


0.14 


0.2 


0.070 


0.22 


0.18 


0.4 


0.085 


0.37 


0.31 



Fig. 7. — Radial structure of strongly magnetized disks in our 
simulations. The upper panel displays disk magnetization, (red 
line with open circles) and dimensionless viscosity generated by 
Maxwell stresses (blue line with open squares). The lower panel 
shows radial dependence of surface density (magenta line with filled 
squared), and mid-plane (red line with open circles) and B r 
(blue line with open squares); B z has the same scaling and magni- 
tude as B r and is not shown here. The green lines show expected 
radial dependence from our analytical model. 

Since the radial dependence of a m is not possible to 
establish from the first principles, we obtain it empir- 
ically. In our numerical experiments a m ex R 1 ^ 2 is 
consistent with the data. This relationship predicts 
that the disk magnetization should decrease with radius, 
f3~ l oc R 2 / 5 ^in agreement with our simulations (upper 
panel in Fig. Y?\. This radial dependence of a m is specific 
to our simulations which we use for consisten cy check. 
Howe ver, a similar dependence was found by |Flock et al.| 
(2011) in the case of weakly-magnetised disks. In a re- 
alistic disk, however, we expect that a m is a function of 
the local /?, which will be subject of subsequent research. 

Using equations above, we derive the radial depen- 
dence of Bp, B r and E. To derive B ri we use the fact 
that in our disks the Maxwell stresses in Eq. Q are dom- 
inated by the mean field, (B r B^) « (B r )(B^). The re- 
sulted radial dependences are 



B Lr 



y _ ( 1600 M 3 C 4 

i 

_ ( 512y/5 M 2 C \ 5 (SIR)* 



(SIR)* „ 7 /oN 
v ; oc R-5, (8) 



ex iJ-s, (9) 



To verify that the mass accretion is physical, we extract 
azimuthally and vertically averaged E, H, M at differ- 
ent radial locations and use Eq.[3] to calculate a acc = 
M/(37rE H 2 Q). If the accretion is driven by magnetohy- 
drodynamical stresses, this value should be comparable 
to azimuthally and vertically averaged sum of Maxwell 
and Reynolds stresses (a). We show results in Tab. [3 
which shows good agreement between measured (a) anc. 
derived viscosity (a acc ) coefficients. This reinforces our 
confidence that the mass accretion is indeed driven by 
magnetohydrodynamical stresses. Finally, using data 
from this table we estimate viscous timescale for a parsec 
size accretion disk to be t v i SC = (R/ i^) 2 a _1 £7 _1 w 10 6 
years, where we use a = 0.2 and H/R = 0.2 at R = 1 pc. 

4. FRAGMENTATION 

The striking result that E and H do not depend on 
the thermal properties of MLAD allows a robust esti- 
mate of its macroscopic gravitational stability^ Its frag- 
mentation boundary is determined by two parameters: 
M and dimensionless mass-to-flux ratio £. The latter 
one allows to compute magnetic flux accretion rate from 
the mass accretion rate, 4> oc M/C, since the magnetic 
field is frozen in the fluid. Using Eq. fidh and Eq . (|8l), 
the Toomre-Q par ameter is ( Toomre||1951| Goldreich~fc| 
Lynden-Bellfl965t 



Q 



Sl 2 H 
ttGE 



6561tt 6 (1 + /3) 3 
64000 G 2 M 2 C 6 



a 2 m (QR) 6 



(11) 



B r = —B^ oc R io. 



(10) 



9 Some of the gas clumps and filaments may form stars even if 
the disk is globally stable. 
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If we write M in terms of Eddington luminosity, Ie 
L/L^dd^ and radiative efficiency, e = L/(Mc 2 ), 



M 



4ttGM l E 



(12) 



where ft es « 0.4 cm 2 /g is electron scattering opacity, and 
c is speed of light, we obtain 



Q 



6561tt 4 k 2 c 2 (1 + /?) 3 < fe_ 



4 5 10 3 



G 2 



(13) 



Using Eq. ( 13 ) we find the fragmentation boundary be- 
yond which Q< 1, 



Rtr: 



2.09 



M 6 a ( 



6 "o.i tp.i 

S10 6 E 



pc, 



(14) 



where we used e = O.leo.i, a m = O.lao.i, C = 10 Go ana 
M = lO 6 MeM0. The radial dependence of enclosed 
mass and H/R within the fragmentation radius is given 

by 



isk 



M 3C 
H _ 3^r 3 / 10 



1Q 9/10 9/10 

« 0.331-^, 

Go 

r 3/10 



2V10C 



0.149 



10 



(15) 
(16) 



where we define r = R/Rfr ag . It is worth noticing, that 
at the fragmentation boundary, Mdi s k and H/R depend 
only on mass-to-flux ratio. 

In future work we will use our MLAD solution to model 
observations of AGN accretion disk. Here, we briefly con- 



lyj 

sider a parsec-sized disk in NGC1068 (|Jaffe et al.p004) 



as an eaxmple. This Seyfert 2 galaxy hosts an ^ 1O 7 M 
SMBH with a disk extending to a distances of ~ 1 pc. 
The observed upper bound for H/R ~ 0.6 and the hydro 
gen column densi ty Nn °° 10 25 cm -2 (Matt et al. 2004 
Kohler fc Li| 2010), and its luminosity is ~ 0.4 L^dd ( |Pier| 
et al.||1994[ ) Here, we assume that this SMBH accretes 
at Eadington rate (Zg ~ 1) with 10% radiative efficiency 
(eo.i ~ 1); we set ao.i = 1- We find that by setting ( = 3, 
we are able to obtain values for disk thickness and col- 
umn density that are consistent with observations. For 
this parameters, the MLAD thickness at the edge of such 
disk is H/R « 0.15. While this is lower than the ob- 
served value, it is plausible that strong magnetic fields 
contribute toward increasing disk thickness. Finally, us- 
ing the enclosed disk mass at this location and assuming 
that the disk consists purely of atomic hydrogen, we es- 
timate Nn ~ 1-3 x 10 25 cm -2 which is consistent with 
the observational data. Furthermore, the fragmentation 
boundary is located at « 50 pc, which indicates that 
a parsec-sized disk is stable to clumping. This MLAD 
model predicts that such disk should have magnetic field 
strength of ~ 100 mG. 

5. CONCLUSIONS 

In this paper we produce from first principles dynam- 
ically stable models of accretion disks in a state of mag- 
netic levitation. We show that such disks are the natu- 
ral outcome of in-fall of a massive magnetized molecular 



cloud onto supermassive black hole. Such magnetically- 
levitating accretion disks (MLADs) enable large accre- 
tion rates due to the large scale-height and a > 0.1. In 
our simulations, the geometry and strength of the large- 
scale magnetic field are stable for at least 0.24 Myr, cor- 
responding to several hundred orbits at the disk inner 
edge. With measured accretion rates of « O.O5M /yr, 
this is more than 10% of the disk lifetime, supporting 
the claim that such magnetic fields structure is possi- 
bly long lived. The viscous time-scale of such magneti- 
cally levitating disk is estimated to be few million years. 
Interestingly, this feature may help so lve a theoretical 



problem that was recently identified by Alexander et al. 
(2011) with respect to the formation of the stellar disc 



in our Galactic Centre. These authors show that if, as 



according to the currently accepted scenario ( Levin fc 
Belo borodov|2QQ3| |Nayakshin et al.|2007| |Bonnell fc Hrce 
2008), the stellar disc formed as a result of fragment a- 



tion of the massive gaseous accretion disc several million 
years ago, then a substantial gaseous remnant of the ac- 
cretion disc should survive to the present epoch, due to 
the expected long viscous time of the standard Shakura- 
Syunyaev thin discs. Such a remnant is not observed. 
On the other hand, the expected short lifetime of the 
MLAD, may solve the problem of the missing remnant 
gas disk. 

A unique property of magnetically-levitating disks is 
that their surface density and scale-height are indepen- 
dent of the disk's thermal structure. This is expected 
because thermal effects are superseded by magnetic prop- 
erties in determining disk structure. Magnetic levitation 
allows the disk to withstand its own self-gravity to large 
distances. A strong dependence of the fragmentation ra- 
dius on the mass-to-flux ratio of the parent cloud per- 
mits a scenario in which a tidal disruption of a magne- 
tized cloud forms a magnetized gas ring. The inner parts 
spreads inwards on a time-scale determined by global 
magnetic stresses which fuel fast accretion onto the cen- 
tral supermassive black hole, while the outer part frag- 
ments into stars. 

A proper understanding of the field confinement re- 
quires both local and global analysis. In accordance 
with JL, the field confinement appears to be a local phe- 
nomenon and its stability is likely to depend on the rel- 
ative strength of vertical and azimuthal magnetic field, 
which itself depends on both kinematics and magneti- 
sation of the infalling matter. However our simulations 
show that there are non-local processes which generate 
a global coherent magnetic field structure. Its topology 
and strength is ver y important for accretion flows near 
black- hole horizons ([B eckwith et al. 2008; Tche khovskoyj 
et al.|2011||Tchekhovskoy & McKmney|2012p . Therefore 
it is also important to understand the long-term evolu- 
tion of the field topology across several decades in the 
disc radius we believe that both local and global simula- 
tions are essential to our understanding of MLADs. 
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Fig. 8. — The circularly polarized Alfven wave after five crossings of computational domain. The panels show the y-component of the 
magnetic field as a function of x-coordinate. The solid line demonstrates exact solution, while the open circles show the result of simulations. 
In the left-most, middle and right-most panels, the wavelength is resolveid with an average of 13, 26 and 52 meshpoints respectively. 
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APPENDIX 



NUMERICAL METHOD 

In this appendix we demonstrate the ability of our numerical scheme to model MHD flows. Our numer ical method 
combines a mo ving-mesh approach (Trease 1988 Springel 2010) with a weighted particle MHD scheme (Gaburov & 
Nitadori 2011). At every time-step a Voronoi mesh is re- built on the set of particles which is used to solve equations 
ol idea l MHD in the same way as in the weight ed particle scheme. Similar approach has also been attempted in 
TESS flDuffell fc MacFadyenJgpII) and AR EPO (jPakmor et al.|[20TT| ) moving-mesh codes. In contrast to these two 
approaches but similarly to Gaburov & lNitadori| ( 2011| ), we add a source term to the induction equation that restores 
Galliean invariance of scheme in the case of V • B ^ U. This proved to be crucial to stabilize the numerical scheme in 
the presence of strongly magnetised super-alfvenic flows. The numerical code is pubcly available pj Here, we present 
validation of our numerical method by means of two test problems: propagation of circularly polarized Alfven wave 
and the linear regime of magneto-rotational instability in a cylindrical disk. 

Circularly polarized Alfven wave 



This problem was first presented by Toth ( 2000[) as an exact non-linear test problem for ideal MHD. Following |T6th 
(2000) we set the following initial conditions! We use a periodic three-dimensional computational domain with the 
total number of particles equal to N to t = N x x N x /2 x 16, where Nx = 16,32,64 and 128. Particles were in itially 
randomly sampled from a uniform distribution and regularized with the Lloyd's algorithm (e.g. Springel| ( 2010| )). The 



initial conditions are p = 1, P gas = 0.1, B x = 1, v x = 1, B y = v y = 0.1 sin(47nc), B z = v z =0.1 cos(47rx), which 
fits two wave-length into the x-direction. With these initial conditions, the wave-length is resolved with an average of 
6, 13, 26 and 52 mesh-points from the lowest to the highest resolution respectively. The apparent discrepancy from the 
expected resolutions of 8, 16, 32 and 64 mesh-cells per wave-length is due to the initial particle distribution is not being 
a simple cubic lattice, but rather a random distribution which was relaxed by the Lloyd's algorithm. This relaxed 
distribution consists of mesh-cells which can be approximated by regular convex polyhedra with large number of faces 
(> 15). The effective size of such mesh-cell can be approximated by the diameter of a sphere having the same volume 
as the cell itself, and this in turn increases the effective size of the mesh-cell by approximately ~ 1.24 compared 

to a simple cubic cell, while keeping the total volume the same. 

In Fig. [8] we compare the simulated results to the analytical solution. It can be seen that lower resolution simulations 
have more dissipation but do not introduce phase error in the solution. The dissipation is not the result of the 
underlying Riemann-solver or a reconstruction method, but rather that of non-linear monotonicity constraints on the 
linear reconstruction model which is required for a stable description of dis continuities. The side effect of these is th e 
constraint is that it forces the scheme to be first order accurate at extrema ( |van Leer) 19 79 Iwasaki fc Inutsuka||2011 ). 

In Fig.[9] we show Li error of the simulations solution as a function of the number ol meshpoints per unit wavelength. 
Here, the L\ error is defined as L\ = 1/N^2 i \ fi — / ex | where sum is carried out over all TV mesh-cells, and \ fi — f ex \ 
is absolute deviation of the value in a cells, /i, from the corresponding exact solution, / ex . The result demonstrates 
that the convergence for this problem is consistent with the second-order scheme. 

Magneto-rotational instability in non- stratified cylindrical disk 

In this problem we study the ability of our code to reproduce analytical growth-rates of axisymmetric magneto- 
rotational instability. Our computational domain consist of the three-dimensional non-stratified cylindrical disk. The 
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Fig. 9. — The L\ error as a function of resolution for circularly polarized Alfven wave. The open circles connected by the red solid line 
show results of simulations, and the blue solid line is expected dependence for the second order scheme oc 0(N~ 2 ). The vertical axis show 
L\ error in the solution as a faction of number of meshpoints, N, per wavelength. 

inner and outer radii of the disk are equal to R = 1 and R = 8, and the thickness of the disk is if = 1. We use periodic 
boundary conditions in z direction, and outflow boundaries at R = 1 and R = 8. The total computational domain is 
a box with size [16.6 x 16.6 x 1]. 

We simulated three models with an average 14, 20 and 28 meshpoints in z-direction. The initial density is set to 
unity, and we used isothermal equations of state with constant sound speed c s = 0.1. The gravitational potential is 
equal to (j) = — 1/R, where R = y 1 x 1 + y 2 , and the initial velocity is equal to the Keplerian velocity. Initially, we set 
a uniform magnetic field in 2 < R < 4 annulus of the disk with such strength that results in fastest growth for n = 2 
mode at R ~ 2. Namely we have, B x = B y = and B z « 0.055/n, where n = 2. In other words, at R ~ 2 the fastest 
growing MRI mode has the wavelength Amri ~ if/2. In this setup, the Amri is resolved with approximated 7, 10 and 
14 meshpoints in low, medium and high-resolution simulations respectively. 



In Fig. 10 



we show time evolution of radial magnetic energy E = B 2 /2 as a function of the number of orbits at R = 2 
for three" different resolutions. All simulations show exponential growth rater after approximately one orbital period, 
and the low resolution simulations shows growth rate w 0.6O, whereas the medium and high resolution simulations 
show growth rate « 0.6517 and w 0.75(7 respectively. 
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Fig. 10. — This figure shows radial magnetic energy in an annulus 2 < R < 2 + 1/42 (vertical axis) as a function of the number of 
local orbits at R = 2 (horizontal axis). The solid red line shows time evolution of radial magnetic energy for low resolution simulation (on 
average 7 meshpoints per Amri), the green dashed and blue dotted lines shows the results for medium (10 meshpoints per Amri) and high 
resolution (14 meshpoints per Amri)- The left and right dotted lines show exponential growth with slopes 0.757 an d 0.657 respectively, 
where 7 = 4-7T. 
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